/*	Interactive Gravitational Simulator
	Project: IGS Core
	File: Simulator.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "Body.h"
#include "Gravity.h"
#include "NBodySystem.h"
#include "Simulator.h"
#include "Tree.h"
#include <cassert>
#include <iostream>
#include <omp.h>

namespace IGS
{

// Constructs a simulator for the NBodySystem
Simulator::Simulator(NBodySystem *system, Integrator *stepper) : System(system),
	Stepper(stepper), _Time(0), _NextTime(0), _EvaluateTime(0), _TotalSteps(0),
	Scheme(StepConstant), Accelerator(NULL)
{
	Initialize();
}

// Constructs a simulator for the NBodySystem
Simulator::Simulator(NBodySystem *system, Integrator *stepper,
	SteppingScheme scheme, Tree *tree) : System(system), Stepper(stepper),
	_Time(0), _NextTime(0), _EvaluateTime(0), _TotalSteps(0), Scheme(scheme),
	Accelerator(tree)
{
	Initialize();
}

// Initializes the Simulator for use
void Simulator::Initialize()
{
	static const int MinAcceleratorSize = 256;
	const int N = System->Size();

	if (N < MinAcceleratorSize && Accelerator != NULL)
	{
		std::cerr << "Accelerator is not optimal for N = " << N;
		std::cerr << ", it will not be used." << std::endl;
		Accelerator = NULL;
	}

	Reset();
	Evaluate();
	EndStep();
}

// Resets the system to a state where gravity can be evaluated
void Simulator::Reset()
{
	NBodySystem &system = *System;
	for (int i = 0; i < system.Size(); i++)
	{
		system[i].Derivative.Zero();
		system[i].Potential = 0;
	}
}

// Evaluate the gravitational force
void Simulator::Evaluate()
{
	NBodySystem &system = *System;

	// Use accelerator for hierarchical evaluation
	if (Accelerator != NULL)
	{
		Accelerator->Evaluate(system);
		_EvaluateTime = Accelerator->EvaluateTime();
	}

	// Do brute-force pairwise evaluation
	else
	{
		const double t1 = omp_get_wtime();
		if (System->Size() < MinParallelSize || omp_get_max_threads() == 1)
			SerialEvaluate(system);
		else
			ParallelEvaluate(system);
		const double t2 = omp_get_wtime();

		_EvaluateTime = (t2 - t1) * 1000;
	}

	// Apply derivatives to the system
	#pragma omp parallel for schedule(dynamic)
	for (int i = 0; i < system.Size(); i++)
		system[i].ApplyDerivatives(system[i].Derivative);
}

// Starts the current integration step
void Simulator::BeginStep(double target)
{
	NBodySystem &system = *System;

	switch (Scheme)
	{
	// Constant time steps: Use the target
	// time for the next time of each body
	case StepConstant:
		_Time = system[0].Time;
		_NextTime = target;

		for (int i = 0; i < system.Size(); i++)
		{
			system[i].NextTime = _NextTime;
			system[i].CollisionTimescaleSq = 1E12;
		}
		break;

	// Shared time steps: Pick minimum of the
	// the collision time scale for the next time
	// of each body
	case StepShared:
		_Time = system[0].Time;
		_NextTime = target;

		for (int i = 0; i < system.Size(); i++)
		{
			double dt = sqrt(system[i].CollisionTimescaleSq);
			
			// Enforce limits
			if (dt < MinStepSize) dt = MinStepSize;
			if (dt > MaxStepSize) dt = MaxStepSize;

			double bodyTime = system[i].Time + dt;
			_NextTime = std::min(_NextTime, bodyTime);
		}

		for (int i = 0; i < system.Size(); i++)
		{
			system[i].NextTime = _NextTime;
			system[i].CollisionTimescaleSq = 1E12;
		}
		break;

	// Individual time steps: Pick time step on a per-body
	// basis. Each body will evolve at it's own rate to make
	// the most efficient use of the available resources.
	case StepIndividual:
		throw "This is not yet implemented";
		break;
	}
}

// Ends the current integration step
void Simulator::EndStep()
{
	NBodySystem &system = *System;
	for (int i = 0; i < system.Size(); i++)
		system[i].EndStep();
}

double Simulator::MinTime()
{
	NBodySystem &system = *System;

	double t = 1E12;
	for (int i = 0; i < system.Size(); i++)
		t = std::min(t, system[i].Time);
	return t;
}

// Evolves the simulation forward in time by dt
void Simulator::Evolve(const double dt)
{
	double const target = _Time + dt;
	int steps = 0;

	while (_Time < target)
	{
		BeginStep(target);

		Stepper->Predict(*System);

		assert(EvaluationCorrectionOrder >= 1);
		for (int i = 0; i < EvaluationCorrectionOrder; i++)
		{
			Reset();
			Evaluate();
			Stepper->Correct(*System);
		}

#ifdef UsePECE
		Reset();
		Evaluate();
#endif

		EndStep();

		System->RequireUpdate();
		System->Update();

		steps++;
		_TotalSteps++;
		_Time = MinTime();
	}

	_Steps = steps;
}

// Gets the currently evolving NBodySystem
std::shared_ptr<NBodySystem> Simulator::GetSystem()
{
	return System;
}

// Gets the accelerator being used
std::shared_ptr<Tree> Simulator::GetAccelerator()
{
	return Accelerator;
}

// Returns the time spent evaluating forces
double Simulator::EvaluateTime() const
{
	return _EvaluateTime;
}

// Returns the current simulation time
double Simulator::Time() const
{
	return _Time;
}

// Returns the number of steps taken on last integration
int Simulator::Steps() const
{
	return _Steps;
}

// Returns the total number of integration steps taken
int Simulator::TotalSteps() const
{
	return _TotalSteps;
}

}